﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Numerics;
using System.Text;
using System.Threading.Tasks;

namespace SharpSoft.Maths
{
    public static class MathHelper
    {
        /// <summary>
        /// 生成指定两个数之间的所有质数（包括这两个数本身）
        /// </summary>
        /// <param name="min"></param>
        /// <param name="max"></param>
        /// <returns></returns>
        public static IEnumerable<BigInteger> GenPrimalities(BigInteger min, BigInteger max, int testTimes = 20)
        {
            for (BigInteger i = min; i < max; i++)
            {
                if (IsPrimality(i, testTimes))
                {
                    yield return i;
                }
            }
            yield break;
        }


        #region 米勒罗宾 - 质数检测

        private static Random random = new Random();

        private static int[] knowPrimalities = null;
        /// <summary>
        /// 判断一个数是否是质数(米勒罗宾算法)
        /// </summary>
        /// <param name="numberToCalculate">要测试的数</param>
        /// <param name="testTimes">测试次数，测试次数越多结果的可信度越高</param>
        /// <returns></returns>
        public static bool IsPrimality(BigInteger numberToCalculate, int testTimes)
        {
            if (numberToCalculate <= 1)
                return false;
            if (knowPrimalities == null)
            { //常见的已知质数，从小到大排列，这些质数直接比较
                knowPrimalities = new int[] { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29 };
                //knowPrimalities = new int[] { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149 };
            }
            if (numberToCalculate <= knowPrimalities[knowPrimalities.Length - 1])
            {//可以直接判断的数
                return knowPrimalities.Contains((int)numberToCalculate);
            }

            if (numberToCalculate % 2 == 0 || numberToCalculate % 3 == 0 || numberToCalculate % 5 == 0 || numberToCalculate % 7 == 0)
            {//使用2357取模直接排除大部分合数
                return false;
            }

            // Find r such that n = 2^d * r + 1  
            // for some r >= 1 
            BigInteger d = numberToCalculate - 1;

            while (d % 2 == 0)
                d /= 2;

            bool result = true;

            int actCount = 0;//记录实际执行检测的次数
            int parallelCount = 25;//每次并行计算的量
            while (actCount < testTimes)
            {
                //并行运行米勒罗宾测试
                Parallel.For(0, Math.Min(parallelCount, testTimes - actCount), (i, x) =>
                              {
                                  actCount++;
                                  if (!MillerRabbinTest(d, numberToCalculate))
                                  {
                                      result = false;
                                      x.Break();
                                  }
                              });
            }


            return result;
        }
        /// <summary>
        /// 米勒罗宾测试
        /// </summary>
        /// <param name="exponent"></param>
        /// <param name="number"></param>
        /// <returns></returns>
        private static bool MillerRabbinTest(BigInteger exponent, BigInteger number)
        {

            // Pick a random number in 2 to n-2
            BigInteger a = 2 + RandomIntegerBelow(BigInteger.Subtract(number, 4));

            // Compute a^d % n 
            BigInteger x = BigInteger.ModPow(a, exponent, number);


            if (x == 1 || x == number - 1)
                return true;

            // Keep squaring x while one of the 
            // following doesn't happen 
            // (i) d does not reach n-1 this means number is not prime
            // (ii) (x^2) % n is not 1 this means number is not prime
            // (iii) (x^2) % n is not n-1 this means number is prime
            while (exponent != number - 1)
            {
                x = (x * x) % number;
                exponent *= 2;

                if (x == 1)
                    return false;
                if (x == number - 1)
                    return true;
            }

            // Return composite 
            return false;
        }

        public static BigInteger RandomIntegerBelow(BigInteger N)
        {
            byte[] bytes = N.ToByteArray();
            BigInteger R;

            do
            {
                random.NextBytes(bytes);
                bytes[bytes.Length - 1] &= (byte)0x7F; //force sign bit to positive
                R = new BigInteger(bytes);
            } while (R >= N);

            return R;
        }
        #endregion

        #region AKS - 质数检测

        /**
       * @brief Step 1 - If n = a^b for integers a > 1 and b > 1, output composite
       */
        static private bool AKS_test1(BigInteger n, ref bool isPrime)
        {

            isPrime = true;

            const double EPSILON = 1E-10;
            for (double b = 2; b < BigInteger.Log(n, 2); ++b)
            {
                var a = Math.Exp(BigInteger.Log(n) / b);
                if (Math.Abs(Math.Floor(a) - a) < EPSILON)
                {
                    isPrime = false;
                    return false;
                }
            }

            return true;
        }

        /**
        * @brief Step 2 - Find the smallest r such that Or(n) > (log2 n)^2
        */
        static private BigInteger AKS_step2(BigInteger n)
        {

            BigInteger maxk = new BigInteger(Math.Floor(Math.Pow(BigInteger.Log(n, 2), 2)));
            BigInteger maxr = new BigInteger(Math.Max(3.0, Math.Ceiling(Math.Pow(BigInteger.Log(n, 2), 5))));
            bool nextR = true;
            BigInteger r = 0;

            for (r = 2; nextR && r < maxr; ++r)
            {
                nextR = false;
                for (BigInteger k = 1; (!nextR) && k < maxk; ++k)
                {
                    BigInteger mod = BigInteger.ModPow(n, k, r);
                    nextR = mod == 1 || mod == 0;
                }
            }
            --r;

            return r;
        }

        /**
        * @brief Step 3 - If 1 < gcd(a,n) < n for some a ≤ r, output composite
        */
        static private bool AKS_test3(BigInteger n, BigInteger r, ref bool isPrime)
        {

            isPrime = true;

            for (BigInteger i = r; i > 1; --i)
            {
                BigInteger gcd = BigInteger.GreatestCommonDivisor(i, n);

                if (1 < gcd && gcd < n)
                {
                    isPrime = false;
                    return false;
                }
            }

            return true;
        }

        /**
        * @brief Step 4 - If n ≤ r, output prime
        */
        static private bool AKS_test4(BigInteger n, BigInteger r, ref bool isPrime)
        {
            isPrime = (n <= r);
            return !isPrime;
        }

        /**
        * @brief Step 5 - check that for every coeficient (ai) in (x-1)^n ai%n == 0
        */
        static private bool AKS_test5(BigInteger n, ref bool isPrime)
        {
            isPrime = true;

            BigInteger current = 1;

            for (BigInteger i = 1; i < n / 2 + 1; ++i)
            {
                current *= n + 1 - i;
                current /= i;

                if (current % n != 0)
                {
                    isPrime = false;
                    return false;
                }
            }

            return true;
        }
        /// <summary>
        /// 通过AKS算法测试一个数是否是质数，这个结果是完全可信的。
        /// </summary>
        /// <param name="candidate"></param>
        /// <returns></returns>
        public static bool IsPrimality(BigInteger candidate)
        {
            if (candidate <= 1)
            {
                return false;
            }
            bool isPrime = true;

            if (!AKS_test1(candidate, ref isPrime)) return isPrime;


            BigInteger r = AKS_step2(candidate);


            if (!AKS_test3(candidate, r, ref isPrime)) return isPrime;


            if (!AKS_test4(candidate, r, ref isPrime)) return isPrime;


            if (!AKS_test5(candidate, ref isPrime)) return isPrime;


            // test5(candidate, ref isPrime);

            // Step 6 -  Output prime
            return isPrime;
        }

        #endregion
    }
}
